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ABSTRACT 

We present high-resolution 3D smoothed particle hydrodynamics simulations of the for- 
mation and evolution of protostellar discs in a turbulent molecular cloud. Using a piecewise 
poly tropic equation of state, we perform two sets of simulations. In both cases we find that iso- 
lated systems undergo a fundamentally different evolution than members of binary or multiple 
systems. When formed, isolated systems must accrete mass and increase their specific angular 
momentum, leading to the formation of massive, extended discs, which undergo strong grav- 
itational instabilities and are susceptible to disc fragmentation. Fragments with initial masses 
of 5.5 Mj up, 7.4Mjup and 12Mjup are produced in our simulations. In binaries and small 
clusters, we observe that due to competition for material from the parent core, members do 
not accrete significant amounts of high specific angular momentum gas relative to isolated 
systems. We find that discs in multiple systems are strongly self-gravitating but that they are 
stable against fragmentation due to disc truncation and mass profile steeping by tides, accre- 
tion of high specific angular momentum gas by other members, and angular momentum being 
redirected into members' orbits. In general, we expect disc fragmentation to be less likely in 
clusters and to be more a feature of isolated systems. 
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1 INTRODUCTION 

The overall paradigm under which star formation occurs remains 
under debate. One hypothesis is that stars are thought to form 
as a collection of low mass fragments in collapsing clumps of 
gas, and then undergo competitive accretion as the fragments 
try to accrete gas from the ir common reservoir jBonnell et al.l 
I2OOII : iBonnell & BatdllOOj) . The main competing hypothesis is 
gravitational collapse , where massive star-forming clumps collapse 
and form multiple cores. Each star forms from the gas that is 
available in its own core, with limited accretion from other ma- 
terial in the parent clump jKrumholz et al.ll2005b . Regardless of 
the star formation paradigm, in the case of solar-type stars, the 
end result is most likely m embership in a binary or mul tiple sys- 
tem touquennov & Mavoi|[l991 : Eggenberger et al.l2004l) . Scenar- 
ios for the formation of such systems abound. Amongst the sim- 

i >lest, most ide alised scenari os are the fissi on of a bar-unstable core 
Durisen et al . 1986; Burke rTet all 1 1997b and the fragmentation 
of centrally condensed, rotating, magnetised cores (lBosslll997ah . 



More complicated hypotheses appeal to the chaos of the cl oud en- 
vironment, such a s core-core collis ions (iTurner et al.lll995h . proto- 
stellar encounters ( Shen et al.l201C ), dynamical capture in unstable 
multiple s ystems (Bate et al .112002), and accretion-triggered frag- 
mentation feonnell" 1994':' Whit worth et al j|l995l : lHennebelle et all 
|l004; Offner et al. 2008; Kratter et al. 20 1^^ 

Star formation scenarios involving the rapid collapse of a pro- 
tostellar core following the loss of support against gravity require 
the formation of a massive ac cre tion disc around th e central object 
fe.g. lVorobvov & Basul2007l and lWalch et al.l2009h . Such discs un- 
dergo a short-lived ^0.1 Myr stage where the disc is massive rel- 
ative to its host (0.1 < Md/M^ < 1) and where gravitational 
instabilities operate to transport mass through the disc onto the cen- 
tral protostar ( Vorobyov & Basu 2007). Massive, accreting, and ex- 
tended protostellar discs have been shown to be susceptible to frag- 
mentation, which could be responsible for a range of phenomena 
such as FU O nonis events and early dust processi ng if the clumps 
are disrupted teolev et al.|[20TQl : lNavakshi n"2010Y and the forma- 
tion o^_substellar companions otherwise dVorobvov & Basui2007l 
an d lBolevll2009h . 
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Given the problem that both the standard core accretion 
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planet formation ti niescale and protost ellar disc lifetimes are typ- 
ically a few Myr jHaisch et all 1200 ih . the idea of creatine si- 
ant planets in a few orbit al times via gravitational instabilities 
was revived JBossI 1997bh an d has s ince been the sub j ect of 



sustained in t erest (IBossI l200l l2008l: iMaver et all |2004 
Pickett et all I2OOOI : IPickett & DurisenI l2007l : iBolev et all 



2007 



2006 



Bolevll2009h . Analy t ical w orks constraining inner disc fragmen- 



tation ("Rafikovl l2005lj2007h . and the short cooling times required 
to form long-lived clumps within ^ 10 AU in simulations of 
protoplanetary discs, al ong with observations of massive plan- 
ets on wide orbits (e.g. iMarois et~aD I2OO8I , I2OIOI) , have lead to 
a shift in focus to outer disc (> 40 AU) fragmentation (e.g . 



Stamatellos et al. 2007; Bolev 20091: IPodson Robinson et al.l2009l : 



Vorobvov & Basu 2010 ). Whether the extended disc models used 



in these studies are similar to discs formed from collapsed molecu- 
lar cloud cores remains to be seen. 

In light of the likely connection between outer disc fragmen- 
tation and early protostellar systems, we examine in this paper 
the formation of discs in detail using 3D SPH simulations of core 
collapse in turbulent molecular clouds. In particular we compare 
the early evolution in mass, surface density, specific angular mo- 
mentum, and disc stability between several different systems un- 
der near-identical conditions. In one set of simulations we compare 
an isolated and a binary system, and in the other we compare an 
isolated system and a small cluster. Previous studies of gravita- 
tional instabilities in binary systems have yielded mixed results, 
with some finding that the perturbing companion hinders frag- 
mentation through disc truncation and tidal heating (lNelsonll2000l : 
IMaver et al ]|2Qq1) and alternatively promotes fragmentation, also 
through tidal perturbations (|Boss 2006). However, the evolved sys- 
tems considered in these studies may not be as susceptible to frag- 
mentation as their protobinary counterparts, which we investigate 
here, due to the enhanced importance of gravitational instabilities 
during protostellar disc formation. 

The paper is laid out as follows: in Section 2 we discuss the 
simulations, initial conditions and the clump identification proce- 
dure. In Section 3 we present the results and analysis. Further dis- 
cussion is provided in Section 4, and our conclusions are laid out 
in Section 5. 



2 THE SIMULATIONS 

All calculations were run with Gasoline (IWadslevetal1l2004l) . a 
parallel implementation of TreeSPH. We employ a fixed number 
N of smoothing neighbors, with the main r uns using = 32. Ar- 
tificia l viscosity is the standard prescription (iGingold & Monaghan 
1983) with q; = 1, /3 = 2, controlled with a Balsara switch jBalsars 
198^ . 



2.1 Initial Conditions 

The cloud in this experiment is spherical and uniform with a ra- 
di us of 0.188 pc, a mass of 50 M©, and a temperature of 10 K, as 
in Bate et al ] |2003l . It was seeded with supersonic turbulent veloc- 
ities and is marginally self-bound. To simulate interstellar turbu- 
lence, the velocity field of the cloud was generated on a grid as 
a divergence-free Gaussian random field with an imposed power 
spectrum P{k) oc . The resulting velocity dispersion cr(/) 
varies as and is consistent with the Larson scaling rela- 

tions (lLarson|[l98lh . The velocities were then interpolated from 
the grid to the particles. Finally, the condition that the cloud be 



marginally self-bound gives a normalization for the global velocity 
dispersion of 1.17 km s~^. We describe the results of two simula- 
tions in this paper, run A and run B. Runs A & B are different ran- 
dom realizations of the same ICs, with the same power spectrum, 
but with different initial phases of the velocity field. As a result the 
first objects that form in run A (which are the only ones that we can 
follow owing to the high computational expense of modeling hy- 
drodynamics at high resolution without employing sinks), are dif- 
ferent from those in run B. In run A we have two sub- simulations, 
one in which a binary system naturally arises, and another in which 
we suppress the production of the secondary. We then compare the 
subsequent evolution of the systems in an otherwise identical core 
and cloud environment. In run B an isolated system and a cluster 
are produced roughly coevally at opposing ends of a filamentary 
structure. We compare the evolution in these systems, having been 
produced in similar environments. 



2.1.1 Run A: The Isolated and Binary Systems 

We describe here how we produced the isolated system in run A. 
The fiducial cloud free fall time is given by tff = ^JZi^jZIGpci — 
1.92 X 10^ yr, with = 1.2x10"^^ gcc~^ being the initial cloud 
density. The evolution of the cloud was followed up to t = 1.2 tff . 

At t = 0.88tff a core, whose collapse we follow in detail, 
has already formed a small central protostar. The protostar accretes 
rapidly via a connecting filament until the filament itself produces a 
fragment. The fragment survives first pericenter with the protostar 
and becomes a companion. 

We produced an isolated system for comparison under near- 
identical conditions by tracing the particles comprising the com- 
panion back in time to t = 0.88tff and stirring them by ran- 
domly exchanging their velocities. This procedure conserves the 
energy and linear momentum but not the angular momentum of the 
stirred particles. There were approximately 5300 particles, totaling 
0.053 M0, that had their velocities perturbed. 



2.2 Thermodynamics 

We model the thermod ynamics of col laps e with a pi ecewise poly- 
tropic equation of state jTohlindl 19821 and lBatel[T998l) : 



K^P p < 10 

K^p^l'^ 10-^3 ^ P, 



(1) 



with /^o chosen so that the sound speed Cs — 1.84 x 10^ cms~^ 
and the subsequent k\ chosen to ensure pressure continuity. Theo- 
retical estimates of temperatures in molecular clouds find that their 
temperatures should range bet ween 5-10 K over densities rang- 
ing from 10~^^-10~^^ Rcc~^ (iLarsonll 19851 : ILow & Lvnden-Belll 
ll976l : lMasunaga & Inutsukall2000h . Observations indicate a some- 
what higher minim um temperature of 8 K, with typical temper- 
atures of 10-13 K (iKirket al.ll2007l) . The equation of state thus 
captures the approximately isothermal behaviour of the cloud in 
the intermediate density regime. At high densities the gas transi- 
tions to being adiabatic with exponent 7 = 7/5. This equation 
of state is a simplification of the internal heating and cooling pro- 
cesses in a molecular cloud and, in the high density regime, will 
minimize the potential for fragmentation. Furthermore, through 
our choice of softening {tg — 2.3 AU) combined with the poly- 
tropic equation of state, we allow the first hydrostatic core to be 
marginally r esolved, while still inhibiting its subsequent dissociat- 
ing collapse (lLarsonll969h . which can be computationally demand- 
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ing to resolve. Three-dimensionsal collapse simulations with more 
sophisticated thermodynamics, and including radiative transfer, but 
of isolated cores with no turbulent cloud environmen t, have been 
perfo r med bv (iWhitehouse & Ba te 2006 ; Bate 2010 : .Tomida et all 
l201Ql : ICommercon et alJl20^ . While proper treatments of both 
thermodynamics and radiative transfer are important, such simu- 
lations are computationally expensive and so far can only be car- 
ried out at lower resolution, and with shorter integration times than 
required in our calculations. We shall explore the problem with 
more detailed thermodynamics and radiative transfer in a future 
publication. We reiterate, though, that the early evolution of pro- 
tostellar discs and protobinaries is still not understood. Before the 
complexities of radiation hydrodynamics can be introduced, the 
birth of these systems must be understood through examination of 
high-resolution hydrodynamics simulations with a reasonable, al- 
beit simplified, EOS, which is the focus of this paper. 
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2.3 Resolution 

The simulations employ 5M particles. The particle masses 
are thus 1.0 x 10~^ in all simulations. The minimum jeans 
mass Mj, defined as = 2.92 cJ/(G'^/Vj^^). is that ob- 
tained at the transition density to the adiabatic equation of state, 
and is M* = 1.7 x 10"^ M© = l.SMjup. Previous work has 
shown that molecular clou ds should have a minimum jeans mass 
JLow & Lvnden-Belll 19761) , with the gas becoming optically thick 
to its own cooling radiation at'-^ 10~^^gcc~^.In our 5 M parti- 
cle runs Mj /mp = 167 and the jeans mass remains adequately 
resolved at all times. A fixed gravitational softening Cg = 2.3 AU 
is used in all simulations. 



2.4 Disc Identification 

In order to compare and separate the systems in the simulations, 
a working definition of what comprises a protostellar system is 
required. We define a protostellar system within this context as 
being a self-bound gas structure with a peak density greater than 
6.7 X 10~^^gcc~^. The density threshold was chosen to be well 
above the critical density in order to identify only bound structures 
that are in the adiabatic regime. To determine whether any cold, 
isothermal gas is bound to a given system, we require that it is 
bound below a threshold binding energy Eb = 2.6 X 10^^ erg. The 
binding energy criterion is imposed in order to exclude particles 
that are only marginally or not uniquely bound to the system, e.g., 
lying in the protostellar envelope or flowing between the primary 
and secondary. The value was chosen empirically to be as small as 

possible while minimising noise in the results. 

To put this into practice we use the SKID dStadefcOOlh group 
finder. SKID works by pushing tracer particles along density gradi- 
ents to find local maxima, and linking them with friends-of-friends. 
Once group assignments are found we then remove unbound parti- 
cles to create self-bound groups, computing: 



Figure 1. The column density and velocity field in a (1000 AU)^ volume, 
containing 0.29 of gas, around the central prestellar object. The time 
ist = 0.88tff , at the onset of core collapse. 



' + U + Et, 



(2) 



as the total energy, with v and Vcm the velocity and center of mass 
velocity, U the gravitational potential, and Eth the thermal energy. 
The center of mass frame is updated throughout the unbinding pro- 
cedure. 
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Figure 2. The protostellar object formed at the onset of its parental core's 
collapse att = 0.88tff. The object is 1.7 kyr old at this stage with a disc 
of 25 AU in radius, and a mass of 0.055 Mq. 



3 Run A: RESULTS 

At t = O.SStff a nascent protostellar system has begun to form 
from the collapse of its enveloping core. In Fig.[T] the surface den- 
sity and velocity field, projected through a (1000 AU)^ box cut 
around the prestellar object, are shown. The box contains 0.29 M© 
of gas. The central prestellar object formed via fragmentation of 
a filament 1000 AU in extent. We see in the integrated transverse 
velocity field that accretion primarily occurs along the filament, as 
there is only a small torus of opening angles around the system from 
which low density material is observed to make a direct approach. 
Most gas collides with the filament and is subsequently funneled 
onto the system. 

In Fig. [21 the surface density of the prestellar system, in pro- 
jection along its axis of rotation, is shown. The system is 0.055 Mq 
and the disc is 25 AU. The filaments feeding the disc at this stage 
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have densities ranging from 7 x 10~^^gcc~^ up to the critical 
adiabatic density 10~^^ gcc~^. The system accretes rapidly until 
t = 0.893tff whereupon a fragment of mass 8 x 10~^ M© forms 
in the low density filamentary material, close to the Jeans mass of 
6.4 X 10~^ M0. The fragment survives first pericenter and becomes 
a companion. The isolated system, in contrast, continued accreting 
material rapidly, building up the disc around the protostellar sys- 
tem. 



3.1 Surface densities and temperatures 

3.1.1 Surface Density Maps 

In this section, we compare the divergent evolution of the iso- 
lated and binary systems. To accommodate the change in dynamical 
times at disc scales we set the units to kyr and count time relative 
to the formation of the system at O.SStff. In Fig. [3] we see the sur- 
face density maps of the isolated and the binary systems in the left 
and right columns respectively, shown at times 4.4 kyr, 5.8 kyr, 
and 7.1 kyr after formation of the primary (or equivalently, the iso- 
lated system). Differences appear early between the two systems. 
In the topmost row, we see that the primary in the binary system 
is similar to the isolated system in morphology and surface den- 
sity except for a tidal disturbance of the disc caused by an initial 
close passage of the secondary. In the central row the differences 
are more pronounced. At 5.8 kyr we see that the isolated system 
has continued to grow in mass and the disc now has a radius of 
50 AU. Between 4.4 kyr and 5.8 kyr approximately 1 orbital time 
has passed in the binary system and the secondary is at its apogee. 
The primary has shrunk in extent as the secondary has accreted 
some the gas in its outer regions. The secondary is also rapidly 
accreting high specific angular momentum gas from the two fila- 
ments feeding the system. In the bottom row at 7.1 kyr we see that 
the isolated system has continued to grow in mass and extent but 
has become gravitationally unstable and an m = 2 spiral arm has 
developed into a material arm. At 7.1 kyr the binary system has 
completed another orbit, both systems have grown in mass, and the 
same pattern of tidal disturbance of the primary is observed. 



3.1.2 Temperature Maps 

In Fig. m mass-weighted temperature maps of the binary and iso- 
lated systems, at times 4.4 kyr, 5.8 kyr and 7.1 kyr, are plotted, as 
in Fig. [3] As we are using a piecewise polytropic equation of state 
(Equation (TJ, some of the gas bound at the outer edge of the sys- 
tem(s) and some interarm gas is still in the isothermal phase. In the 
isolated system, increasing amounts of mass end up at large radii 
and at low densities, due to its higher specific angular momentum j 
(see section [T2I for elaboration on these points). In the binary sys- 
tem however, pictured in the right column of Fig. U almost all of 
the mass identified as bound to each system lies at high densities 
(and is thus adiabatic) and within their respective tidal radii. 



3.2 Disc Evolution 

3.2.1 Disc Mass 

In Fig. [5] the accretion histories of the isolated system, primary, and 
secondary are plotted. The mass at a given time is determined by 
the procedure outlined in 12.41 and is thus the combined protostel- 
lar and disc mass. The isolated system accretes gas steadily and at a 
high rate throughout the simulation, with typical accretion rates of 



O 
C\J 



if) 

CO o 
E d 



isolated 




— — primary 

primary + secondary 


















^ ^ . 






- 

^0* 





3 4 5 6 7 
time (kyr) 



Figure 5. Total mass accreted (disc + protostar) versus time for the isolated, 
primary, and secondary systems, and the net mass in the binary system. 
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Figure 6. The ratio of disc mass to accreted stellar mass versus time, for 
the isolated, primary systems. 



3 — 5 X 10~^ M0 yr~^. After the initial collapse, the central pro- 
tostar, defined as the mass contained within r < 5 AU 2e, has 
an accretion rate of 5 — 6 x lO~^M0yr~^. Taking the asymp- 
totic rate from ID pr otostellar collapse theory, M — moCg/G 
(IStahler & Pallall2005h . we find that M = 2 x 10"^ M© yr-\ and 
so TTio, a constant of order unity, is approximately 2.5 at the pro- 
tostar. This is not a high accretion rate, as ID collapse simulations 
asymptote to mo = 2 (IStahler & Pallall2005h . The accretion rate 
onto the entire system however, is high, with vfio ^ 10, and is per- 
haps due to embedding within a filament. Both accretion rates are 




Figure 3. The surface densities of the isolated and binary systems, respectively in the left and right columns, and shown at times 4.4 kyr, 5.8 kyr, and 7.1 kyr 
from the top frow to bottom, after formation. 
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Figure 7. Evolution of the surface density profile of the isolated system, a 
tends to increase at all radii with time. 



consi stent with other studies, such as those found in IWalch et al J 
I2OO9I . 

As expected, the net mass accreted by the binary system is 
close to that of the isolated one. The difference between the two 
is because there is some mass present in the binary that can not 
be uniquely assigned either to the primary or secondary. The bi- 
nary system is plotted starting at 3.8 kyr, instead of its formation 
time at 3330 yr, because during the first pericenter the identifica- 
tion procedure has difficulties separating the two objects uniquely. 
In the binary system the secondary, due to its orbit, preferentially 
accretes high specific angular momentum gas from the filaments 
and accretes mass at a higher rate than the primary, and at later 
times the mass ratio tends towar d unity. This behaviou r has been 
noted before in the literature fe.g. lBate & Bonnelllll997h . 

In Fig. [6] the ratio of stellar mass (the mass within 5 AU, cho- 
sen to correspond with the limits of our softening) to disc mass for 
the isolated and primary systems are plotted. The isolated system 
continues gaining disc mass faster than it can be transported to the 
star, while the primary evolves steadily towards stability. As sys- 
tems normally become moderately self-gravitating at a mass ratio 
of 0.1, we see that these newly formed protostellar systems are 
in a regime dominated by gravitational instability. 
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Figure 8. Evolution of the surface density profile of the primary in the 
binary system, a decreases steadily at outer radii and increases at inner 
radii due to tidal interactions. 
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3.2.2 Cumulative Mass Fraction and Mass Profiles 

The discs in both the isolated and binary systems are all massive 
relative to the central protostar and exhibit strong spiral arms as 
seen in Figs. [3] We expect that tidal torques will, in addition, play 
a role in triggering mass transport within th e discs in the binary 
as noted for example in jMaver et al.ll2005h . In Figs. |7]&[8]are 
plotted the surface densities averaged in annuli of the isolated disc 
and of the primary of the binary system. We see from Fig. [7] that 
in the isolated system the trend as time progresses is toward in- 
creasing amounts of mass at all radii, while from Fig. [9] it is ev- 
ident that the relative amount of mass at large radii is increasing, 
illustrated by the outward movement of the half-mass radius. The 
surface density profile of the primary in the binary system (see 



Figure 9. The mass profile evolution of the isolated system. As time pro- 
gresses, increasing amounts of mass relative to the total lies at large radii. 
We illustrate this by plotting the evolution of the half -mass radius in black. 



Fig. [S]) shows increasing surface density within 20 AU and clearly 
decreasing surface density outside of 30 AU throughout the sim- 
ulation. This is reflected in the cumulative mass profile in fig. [TOl 
as a decreasing overall fraction of mass found at larger radii, and 
the clear inward movement of the half-mass radius (while the total 
mass slowly increases). The source of this difference is likely in- 
creased mass transport from gravitationally amplified tidal pertur- 
bations. To illustrate the role of tidal effects we model the primary 
and secondary as point masses and compute the Jacobi radius Vj at 
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Figure 10. The mass profile evolution of primary. In contrast to Fig. |9] 
we see the relative fraction of mass at large radii decreasing with time. We 
illustrate this by plotting the evolution of the half-mass radius in black. 



20 40 60 80 100 
r (AU) 

Figure 11. Azimuthally mass-averaged temperatures in the isolated system. 
At later times the high j gas has formed an extended disc which includes 
cold interarm gas around the primary, leading to a break in the temperature 
profile. 



5.8 kyr given the nominal values M = 0.125 M©, m = 0.05 M© 
and Rq — 75 AU (see Figs.[5]&[3]) for the primary mass, secondary 
mass and separation, respectively Given these parameters we com- 
pute Tj = 0.41i?o = 31 AU, thus we expect tidal effects to, at 
the very least, limit the disc of the primary to 44 AU. Looking at 
the cumulative mass profile, we see that at our disc-finding proce- 
dure is roughly consistent with such a simple model, with very little 
bound mass found outside of 30 AU. 

3.2.3 Temperature Profiles 

In Figs. [TT]&[T2l are plotted the temperature profiles in the isolated 
and binary systems, respectively At times 4.4 kyr and 5.8 kyr, the 
temperature declines steadily from 250 K to 10 K at 40 AU and 
^ 60 AU respectively, whereupon it remains at the temperature 
floor for larger radii, however there is little mass at these radii (see 
Fig-HJ- At time 7.1 kyr, the disc becomes more extended, and the 
temperature profile flattens at large radii. The disk contains a mix- 
ture of isothermal and adiabatic gas, with the spiral arms containing 
most of the adiabatic component. In the binary system (Fig.[I2]), the 
temperature profiles of the primary lack the low temperature, low 
density component seen in the isolated system due to truncation of 
the disc at the tidal radius, and do not evolve much after 5.8 kyr. 

3. 2. 4 Specific Angular Momentum 

In Fig. [13] are plotted the combined disc and protostellar specific 
angular momenta for the isolated system, primary and secondary 
Each measure is computed in the centre-of-mass frame of the given 
systems and is thus a measure of the spin (i.e. disc) component 
of the specific angular momentum. The isolated system rapidly in- 
creases in specific angular momentum, and after 7.1 kyr, becomes 
unstable to fragmentation. In the binary system the primary has a 
small spike in its specific angular momentum j at 4.2 kyr due to 
weakly bound material appearing to be bound to the primary while 
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Figure 12. Azimuthally mass-averaged temperatures in the primary of the 
binary system. The disc remains limited by tidal forces, and almost all 
bound gas is adiabatic. 



the secondary reaches pericenter. This is also seen in the mass plot 
at the same time. The primary exhibits very little evolution in j 
thereafter, while from Fig.|5]we see that it continues to accrete mass 
in the interim; it accretes at constant j. In contrast the secondary in- 
creases steadily in j along with M. 

Part of the reason that the binary increases more slowly in 
specific angular momentum than the isolated system, is that it is 
able to store momentum in the binary orbit. To measure this effect. 
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Figure 13. Evolution of the specific angular momentum j of the isolated 
and binary systems, j grows rapidly for the isolated system and leads to a 
more rapid buildup of an extended disc, whereas tidal effects and lower j 
leads to smaller discs in the binary system. 
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Figure 14. Azimuthally averaged Toomre Q for the isolated system at times 
4.4 kyr, 5.8 kyr and 7.1 kyr. 



we use the relationship 



LO 



(3) 



where L is the angular momentum, are the positions of each 
particle, Pi are the particle momenta, R is the barycentre, P the 
velocity of the barycentre, Vi and positions and momenta rela- 
tive to the barycentric frame, and decompose the binary into spin 
and centre of mass angular momenta. We further get the spe- 
cific angular momenta using the appropriate mass-weightings along 
with Equation [3] At t = 7.1 kyr, we find that in the barycen- 
tric frame, the isolated system has a specific angular momentum 
of 7.0 X 10^^ cm^ s~^. In the binary system, the binary orbit it- 
self has 6.1 X 10^^ cm^ s~\ the primary has 3.0 x 10^^ cm^ 
and the secondary 1.3 x 10^^ cm^ (also visible from Fig.[T3l). 
The specific angular momentum of the binary orbit dominates that 
found in the other components. 

We quantify alignment of the specific angular momentum 
of the disc with that of the environment by computing the angle 
Q — cos ~^(jd • je), with jd being the unit vector of the specific 
angular momentum of the disc and je being that of the material in 
the environment in a box of dimensions 2000"^ AU^ surrounding 
the system. At the time of formation of the system (Figs. [T]&[2l) 
we find Q — 36°. In the isolated system we observe angles of 13°, 
7.7°, 12° at times 4.4 kyr, 5.8 kyr and 7.1 kyr, respectively. In 
the binary system the jd of the primary is 18°, 7.3°, 5.3° at times 
4.4 kyr, 5.8 kyr and 7.1 kyr. The vectors tend to be oriented per- 
pendicular to the major axes of the filaments feeding the disc. Such 
relatively small misalignments are consistent with other works sug- 
gesting filaments tumbling about one of t heir shorter axes a s the 
origin of protostellar angular momentum dWalch etal.ll2009l) and 
with observations of protostellar systems on these length scales, 
which find t hat most gas is ex tended in the direction perpendicular 
to outflows (iTobin et al.ll201Ql) . 
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Figure 15. Azimuthally averaged Toomre Q for the binary system at times 
4.4 kyr, 5.8 kyr and 7.1 kyr. The binary system remains stable for as long 
as the simulation was run (20 kyr). 

3.3 Stability and Fragmentation 

Throughout the simulation(s), the discs of both the isolated and bi- 
nary systems are massive and gravitationally unstable. To compare 
stability we compute the Toomre Q — CsI^/tiGct parameter locally 
using mass- weighted projections, through a 200 x 200 x 200 AU^ 
box, followed by azimuthal averages, of the disc. The results are 
shown in Figs. [14] & [15] Although the Q curve for the isolated 
system drops below unity in some regions at 4.4 kyr it does not 
fragment right away. As the envelope is accreted, the disc evolves 
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Figure 16. Toomre Q paramater of the isolated system averaged in an 
annulus from 60-70 AU in times preceeding clump formation. A clump 
later forms at 9.2 kyr at 65 AU. 
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Figure 17. A surface density map of the isolated system at t = 11 kyr, 
approximately one orbital time after the second fragment formed. 



to a state where Q lies between one and two out to 90 AU. The 
disk remains in the unstable regime, but Q varies substantially as 
the system accretes. In contrast, the primary of the binary system 
has a minimum Q of 1.3 at later times and, except initially, rises 
above two at radii greater than 25 AU. The binary system is thus 
significantly more stable than the isolated one, and this is borne out 
by the fact that it does not fragment at later times (up to 20 kyr). 

At 9.2 kyr we find that the isolated system becomes unsta- 
ble to fragmentation, forming two clumps, one at r = 65 AU, 
followed by another at r = 100 AU in the disc. The first clump 
(CI) formed in one of the material arms of the disc from adia- 
batic gas, with an initial mass of 5.5 Mj up and specific angular mo- 
mentum of 1.8 X 10^^ cm^ s~^. The second clump (C2) formed at 



10 kyr, had a mass of 7.4 Mjup and specific angular momentum of 
3.2 X 10^^ cm^ s~^. We seek to understand formation of the first 
clump from adiabatic gas by looking at the behavior of Q in the 
run-up to fragmentation. In Fig. [161 is plotted the Q parameter av- 
eraged in an annulus ranging from 60-70 AU. Q starts off initially 
very high but drops as continuing accretion builds up the disc in 
this region. Both clumps are depicted in Fig.[T7]at 11 kyr. At this 
time, clump CI has made approximately 1.5 orbits and has a mass 
of 39 Mjup, and clump C2 one half-orbit, has grown to 14 Mjup. 
Both clumps accrete at a rate of roughly 10~^ M© yr~^. 



4 Run B: RESULTS 

We describe here the results from a complimentary simulation, run 
B. At approximately t — 1.1 tff fragments form at either end of 
a long (7400 AU) filament, within 2 kyr of each other. In the 
following discussion and figures, we set t = 1.1 tff as being t = 0, 
as it is at this time when the first system, the isolated system, forms. 

The upper left image in Fig. [18] is a surface density plot of 
the filament, projected through a (10^ AU)^ box. At the top end 
of the filament lies the isolated system and the cluster is located at 
the bottom. The upper right image in Fig.[T8]is the surface density 
of the cluster at the end of simulation. Only three of four members 
are shown. The two systems in the upper left corner of the image 
have formed a binary, in the sense that they complete orbits around 
their centre of mass on a time scale shorter than the cluster crossing 
time. For a more complete discussion of the cluster evolution, see 

El 

In the second and third rows of Fig. [18] are the surface den- 
sities and mass-weighted temperature maps of the isolated system 
(lying at the top of the filament), at 4.3 kyr, and 8.7 kyr. In the 
isolated system we see the disc assemble rapidly from infalling fil- 
amentary gas, going through stages of increasing gravitational in- 
stability, until, as shown in the bottom row of Fig. [18] gas trailing at 
the end of one of the material arms in the disc fragments. As seen in 
the temperature maps, the gas in the fragment was in the adiabatic 
regime. The initial clump mass is 12 Mjup. 



4.1 Cluster Evolution 

In Fig. 14. II are plotted the distances of cluster members to the ori- 
gin of a fiducial cluster frame vs. time. Since the cluster is a rapidly 
evolving object at this stage, we choose the fiducial frame as being 
the one which moves in a straight line from the original position 
of the first clump to its final position. We can see in Fig. 14. II that 
members are typically 100 AU from the origin in this frame. 
Members spend most of their time in binary associations, whereby 
it is meant that they orbit around their common center of mass on 
timescales shorter than the cluster crossing time. Two such binaries 
can be seen at points A and B in Fig. 14. II through the mutual oscil- 
lation of their distances in the cluster frame. At point C in the figure, 
the members of the cluster converge briefly and their mutual sepa- 
rations become comparable to their binary separations. Between C 
and D, a merger occurs, and a new binary forms, with the remaining 
member being sent on a larger orbit. 
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Figure 18. First row, Left: Surface density map of filament formed in run B. At the top and bottom of the filament are the isolated system and cluster, 
respectively. Right: The cluster as seen at the end of the simulation. One member, 10^ AU away, is not shown. Second and third rows: The surface densities 
and mass-weighted temperature maps of the isolated system, respectively in the left and right columns, shown at times 4.3 kyr, and 8.7 kyr, after formation. 
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Figure 19. Distance from the origin of a fiducial cluster frame of reference 
vs. time, for various members of the cluster. The frame of reference is based 
on the trajectory of the oldest member, hence the separation approaches zero 
at the beginning and end of the data. Several objects pair up to form short 
lived 'binaries' while in the cluster. Two such binaries can be seen at points 
A and B in Fig. 14. II through the mutual oscillation of their distances in the 
cluster frame. At point C, all extant members of the cluster briefly converge. 
Between C and D, a merger occurs, leading to the formation of a new binary 
and the ejection of the third member. 
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Figure 20. Total mass accreted (disc + protostar) versus time for the iso- 
lated system, and the cluster. A merger occurs between two prestellar sys- 
tems in the cluster at 6.3 kyr. 
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Figure 21. Evolution of the specific angular momentum j of the isolated 
system and the cluster, j grows rapidly for the isolated system and leads to 
a more rapid buildup of an extended disc, whereas various tidal effects and 
close encounters lower j and leads to smaller discs in the binary system. 

4.2 Disc Evolution 

4.2.1 Disc Mass 

In Fig. |20l the accretion histories of the isolated system, and those 
of all objects formed in the cluster are plotted. The isolated system 
accretes gas at a high rate, at an average of 8 x 10~^ M© yr~^. The 
accretion rates observed are quite comparable to those for the iso- 
lated system from run A, and thus they are well above the expected 
asymptotic rates (mo > 10). 

For the cluster, the average accretion rate is 3 x 
lO~^M0yr~^, with little deviation between objects. The evolu- 
tion of the cluster is complex. The pattern of fragmentation of core 
gas never yielded objects which were isolated for very long, quickly 
leading to the production of a couple of short-lived binaries with 
similar evolution (similar accretion rates and generally limited disc 
size). A major merger occurs at 6.3 kyr; the result of a head-on 
collision between two clumps. Given the chaos of the cluster, the 
accretion rates as seen by the slopes in Fig.|20]are remarkably uni- 
form. 



4. 2. 2 Specific Angular Momentum 

In Fig. [21] are plotted the combined disc and protostellar specific 
angular momenta for the isolated system, and for the objects in the 
cluster. As in run A, the isolated system rapidly increases in specific 
angular momentum. Upon reaching ^ 7 x 10^^ cm^ s~^, material 
arms form in the isolated system, moving angular momentum out 
of the system, and allowing angular momentum to remain constant 
while accretion continues (compare to Fig. |20]). As the arms be- 
come stronger, the determination of j becomes more difficult, be- 
cause gas is ejected with increasing violence from the ends of the 
material arms back into the envelope. About 1 kyr after the onset 
of this behaviour, the arm fragments (shown in Fig.fTS]). 

The evolution of the cluster angular momenta shows a rather 
remarkable contrast. Following Fig. |2T] once more, we see that the 
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first object to form in the cluster evolves as one might expect, with 
both mass and specific angular momentum increasing. Once other 
fragments form in the cluster, the specific angular momentum of 
the system drops until hitting a maximum value which appears to 
characterise the entire cluster. The pattern of fragmentation in the 
core never left any object in the cluster truly isolated for very long, 
with some short-lived binaries forming, sometimes ending in merg- 
ers, but generally lasting long enough to tidally limit and alter the 
structure of the discs of cluster members. Referring once more to 
the mass evolution plots in Fig|20l we see that accretion continues, 
somewhat abated by the competition amongst cluster members for 
inf ailing gas. We see then a stronger, clearer, reiteration of the re- 
sult from run A, namely that upon hitting a threshold, accretion 
continues but at constant specific angular momentum. 

The major merger occuring at 6.3 kyr, indicated by the large 
filled circles in Fig.|2T]has an interesting effect on the angular mo- 
menta of the merging systems. Their momenta both dip quite sud- 
denly at the time of collision. This is due to the ejection of plumes 
of high-j gas at the moment of merger. The gas is then quickly 
accreted onto the new system, bringing it back up to the angular 
momentum ceiling for the cluster. 



5 DISCUSSION 
5.1 Fragmentation 

5.1.1 Accretion and specific angular momentum 

In the simulations performed the isolated systems become prone to 
fragmentation via the development of massive material arms, even 
though the gas is adiabatic. Instabilities in the isolated systems are 
driven by steady accretion of gas with increasing specific angular 
momentum. The disc, sufficiently massive to generate spiral arms, 
accumulates a reservoir of gas lying at large radii (i.e. with high 
specific angular momentum). This reservoir is then swept up by 
the spiral arms, which when sufficiently massive, become prone to 
fragmentation; such instabili ties have been documented in similar 
contexts in the literatur e re.g. lBonnelllll994|:IWhitworth et al.ll995l : 
iHennebelle etaD l2004^. At the same timeJOffne r et all bOOSh ob- 
serve that in simulations of freely decaying turbulence the resulting 
protostellar systems are more prone to fragmentation, and attribute 
it to high rates of accretion. Seeing as the isolated systems in the 
present study have accretion rates several times that of the nomi- 
nal ID asym ptotic rate, perhaps the tendency toward fragmentation 
observed bv lOffner et al.., ( 2008) is due to a combination of higher 
accretion rates and accretion of gas with high specific momentum. 
Furthermore, it could be that the continued driving of turbulence 
is reducing the specific angular momentum of accreted gas, or act- 
ing as a viscosity and depleting the reservoirs of gas which would 
otherwise be swept up in dense material arms. 

The binary system from run A and the systems in the cluster 
from run B remain stable throughout the simulation. Although the 
discs are sufficiently massive to be self-gravitating, and the systems 
are accreting, a number of effects work against outer disc fragmen- 
tation. The tidal limiting of the discs, as discussed in section [T2l2l 
is a major barrier to fragmentation. In run A, although the binary 
separation evolves (increases) throughout the simulation, by the 
end of our run the separation is ^ 75 AU, a nd given the tendency o f 
the binary mass ratio to approach unity fe.g. lBate & Bonnelll 19971) . 
there is very little chance for outer disc fragmentation to occur. In 
run B the cluster appeared to have an effective ceiling in the spe- 
cific angular momentum for member systems. Longer simulations 



looking at the long term evolution of multiple systems and clus- 
ters should be undertaken to see if they separate sufficiently to al- 
low extended discs to grow. Even if at some later point they do 
grow, fragmentation is less likely because most mass will likely 
have made it onto the star and a high disc/star mass ratio will be 
harder to achieve, also, accretion rates decrease as time progresses. 
In run A, the tendency for the secondary to accrete high specific 
angular momen tum gas as seen in thi s paper and in other studies 
of binaries fe.g. lBate & Bonnellll 19971) . acts strongly to deplete the 
reservoir of gas in the outer part of the system and shuts out the pos- 
sibility for the growth of the same Toomre instability which lead to 
fragmentation of the isolated system. Furthermore, the enhanced 
inward transport of mass due to tidal torques shown in l3.2.2l acts to 
stabilise the disc against fragnientat ion, an effect that has also been 
observed by (e.g. lMaver et al.ll2005b . 



5.1.2 Accretion and gas temperatures 

In the following we discuss analyse run A in more detail. In the iso- 
lated system two clumps form by the end of the simulation. They 
form in the outer reaches of the disc where temperatures were low, 
and where some interarm gas is still isothermal at 10 K, although 
clump CI forms in a material arm from adiabatic gas with a tem- 
perature of 30 K. Recent simulations of star formation including 
a flux-limited diffusion treatment for radiative transfer and some 
modelling of protostellar accretion emphasise the importance of ra- 
diative feedback on the envelope. To esti mate the effect we comp ute 



the accretion luminosity of the protostar 



gate the eiiect we comp u 
dStahler & Pallall2005f) : 
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with =0.11 M0 being the protostellar mass, which we com- 
pute by taking all mass within r <= 5AU ^ 2e; M = 5 x 
10~^ M© yr~^ is the (time- averaged) accretion rate onto the pro- 
tostar, and R^ the protostellar radius that, since it is sub-grid, we as- 
sume a fiducial value of 5i?0 . From this we get a typical accretion 
luminosity of Lace — 3.4 L©. Clump CI forms at Vc — 65 AU. 
To try to get an upper limit on the temperature at that radius, 
we assume that the disc presents an absorbing surface of height 
h — 10 AU spanning rc db 5 AU, and that all of the accretion lumi- 
nosity crossing this surface deposits all of its energy in this region, 
yielding a heating rate of Lacch/2rc — 10^^ ergs~^. We assume 
iD'Alessio et al.l (l200lh opacities and a 1 /im grain size, and a back- 
ground irradiation temperature of 30 K, and using the existing den- 
sity field of the gas within rc db 5 AU at the time of fragmentation, 
we calculate the temperature at which the cooling rate balances the 
heating rate of the gas to be 50 K. We wish to emphasise that the 
assumptions made here are extremely conservative with the aim of 
bracketing the upper end of conceivable outer disc temperatures. 
Given the weak dependence of Q on T, it does not appear that the 
accretion luminosity will be able to affect the disc stability at these 
large radii. 



5.1.3 Clump masses and accretion rates 

Again, we limit the discussion here to run A for brevity. The iso- 
lated system forms two clumps of masses 5.5 Mjup and 7.4 Mjup. 
What kind of clump masses should one expect in the outer disc? 
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Because the disc is in a disordered state at the time of clump forma- 
tion we do not compare the clump masses to those computed from 
axisymmetric models, but instead we compute the local Jeans mass 
of the clump-forming gas, which for clump CI was 4.7Mjup, and 
that of clump C2 was 4.5 Mjup. The particle mass in the simula- 
tions was 10~^ M0, hence in the clumps the Jeans mass is resolved 
with ^ 450 particles, or ^ 14 smoothing kernels. The fragmenta- 
tion phenomena are t hus well-reso lved according to the criteria set 
forth in the literature (iBate & Burkert l997 : Nelson 2006) . In addi- 
tion, the initial angular momentum radii of the clumps are 0.31 AU 
and 0.73 AU, respectively, and should they undergo a second col- 
lapse in ^ 10^ years due to H2 dissociation, any accret ion onto the 
final body will take place through an accretion disc (iBolev et aP 
We observe final clump masses of 39 Mjup and 14 Mjup 
for CI and C2, respectively, however whether or not dissociative 
collapse will occur before these masses are achieved is highly sen- 
sitive to assumptions about the dust. We emphasise here that the 
connection between our final obtained clumps and the substellar 
companions (the planets or brown dwarfs) produced from them is 
uncertain. A dissociating collapse of the clumps and the establish- 
ment of an accretion disc is a certainty, however other dynamical 
instabilities may form, such as a bar instability, and can play a role 
in redistributing mass and angular momentum. Because the simula- 
tion lacks the necessary physics and resolution to follow the inter- 
nal evolution of the clumps, we refrain from speculating on the final 
object masses. Nevertheless, we can estimate accretion luminosi- 
ties for our final clumps using our clump parameters. Both clumps 
CI and C2 had accretion rates of'-^ lO~^M0yr~^. Computing 
the accretion luminosity using Equation |4l inputting the measured 
accretion rate, clump masses of 5.5-39 M©, and assuming a fidu- 
cial final size of 1 Rjup = 7 x 10^ cm, gives accretion luminosities 
ranging between 16-110 L©. Such relatively high luminosities are 
achieved at this accretion rate because the radius of Jupiter is ^ 50 
times smaller than our fiducial one for the central protostar (5 R©). 
With such extremely high luminosities we can see that mass ac- 
cretion simply cannot continue unabated without feedback playing 
some role. Our final clump masses and luminosities thus fall on the 
upper end of the range of expected values. 



6 CONCLUSIONS 

The primary motivation of this paper was to compare the evolution 
of an isolated and binary/multiple systems during the early stages 
of prestellar development within the same environment. Of partic- 
ular interest was the potential for fragmentation of such early, ac- 
creting and self-gravitating discs. The following are our findings: 

• The initial collapse occurs within filaments that tend to col- 
limate infalling gas onto the central prestellar objects, whose spe- 
cific angular momenta tend to remain aligned to that of the envi- 
ronment (and perpendicular to the embedding filament) on scales 
of - 1000 AU. 

• During formation, isolated systems must accrete mass and in- 
crease their specific angular momentum, leading to the formation 
of massive, extended discs, which undergo strong gravitational in- 
stabilities and are susceptible to fragmentation 

• In systems starting out as binaries and multiples, an effective 
specific angular momentum ceiling exists, limiting the maximum j 
of the systems but not their mass, making disc fragmentation un- 
likely. The ceiling is the product of tidal interactions which strip 
cold gas from members, and concentrate their mass profiles, and 
the redirection of angular momentum into orbits. 
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